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Abstract 

A  comprehensive  non-isothermal,  3D  computational  model  for  proton  exchange  membrane  (PEM)  fuel  cells  has  been  developed,  and 
implemented  into  a  computational  fluid  dynamic  (CFD)  code.  The  model  allows  parallel  computing,  thus  making  it  practical  to  perform 
well-resolved  simulations  for  large  computational  domains.  The  model  accounts  for  convective  and  diffusive  transport  and  allows  prediction 
of  the  concentration  of  species.  Distributed  heat  generation  associated  with  the  electrochemical  reaction  in  the  cathode  and  anode  is  included 
in  the  model.  The  model  solves  for  the  electric  and  ionic  potentials  in  the  electrodes  and  membrane,  and  the  local  activation  overpotential 
distribution  is  resolved,  rather  than  assumed  uniform,  making  it  possible  to  predict  the  local  current  density  distribution  more  accurately. 

Maximum  current  densities  are  predicted  under  the  land  areas  as  a  result  of  the  dominant  influence  of  ohmic  losses  over  concentration 
losses  on  the  activity  at  the  catalyst  layer.  A  parametric  analysis  shows  that  substantially  different  spatial  distributions  can  be  obtained  by 
varying  the  asymmetry  parameter  with  no  noticeable  change  in  the  global  current  density  and  polarization  curve.  Changing  the  conductivity 
radically  alters  the  current  distribution  by  changing  the  relative  influence  of  ohmic  to  activation  overpotentials. 

©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Recent  years  have  seen  significant  increase  in  power  den¬ 
sities,  reliability  and  overall  performance  of  proton  exchange 
membrane  (PEM)  fuel  cells,  but  the  underlying  physics  of  the 
transport  processes  in  a  fuel  cell —  which  involve  coupled 
fluid  flow,  heat  and  mass  transport  with  electrochemistry — 
remain  poorly  understood.  The  development  of  physically 
representative  models  that  allow  reliable  simulation  of  the 
processes  under  realistic  conditions  is  essential  to  the  devel¬ 
opment  and  optimization  of  fuel  cells,  the  introduction  of 
cheaper  materials  and  fabrication  techniques,  and  the  design 
and  development  of  novel  architectures. 
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The  difficult  experimental  environment  of  fuel  cell 
systems  has  stimulated  efforts  to  develop  model  that  could 
simulate  and  predict  multi-dimensional  coupled  transport 
of  reactants,  heat  and  charged  species  using  computational 
fluid  dynamic  (CFD)  methods.  The  first  applications  of 
CFD  methods  to  fuel  cells  focused  on  two-dimensional 
models  ([1-4]).  More  recently,  CFD  and  improved  transport 
models  have  allowed  the  development  of  increasingly  more 
realistic  computational  models,  accounting  for  fluid,  thermal 
and  electrochemical  transport,  complex  three-dimensional 
geometries  including  flow  and  cooling  channels,  and 
two-phase  transport  (e.g.  [5-8]). 

Simulations  of  3D  PEMFC  geometries  have  required  some 
simplifications  in  order  to  reduce  computational  require¬ 
ments.  In  particular,  it  is  only  very  recently  that  CFD-based 
models  have  started  to  include  resolved  catalyst  layers  and  to 
account  for  non-uniform  distributions  of  overpotentials  at  the 
electrodes  [9-11].  A  common  issue  to  many  computational 
models  is  the  uncertainties  associated  with  the  specification 
of  various  parameters  that  impact  the  transport  processes.  As- 
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sessing  the  sensitivity  of  the  flow,  thermal  or  electrochemical 
response  to  these  uncertainties  is  essential. 

In  this  paper  a  single-phase,  three-dimensional,  non- 
isothermal  model  is  presented.  This  model  is  implemented 
into  the  commercial  CFD  code  Fluent  6.1,  with  custom- 
developed  user-subroutines  that  take  account  of  the  physico¬ 
chemical  processes  associated  with  PEM  fuel  cells.  The 
model  has  the  capability  of  resolving  the  catalyst  layer,  and, 
in  contrast  with  most  published  models,  this  model  accounts 
for  a  distributed  overpotential  at  the  cathodic  catalyst  layer, 
and  heat  sources  at  each  electrode,  rather  than  lumping  the 
heat  sources/sinks  at  the  cathode  as  is  common  practice.  This 
model  also  takes  into  account  convection  and  diffusion  of 
different  species  in  the  channels  as  well  as  in  the  porous  gas 
diffusion  layer,  heat  transfer  in  the  solids  as  well  as  in  the 
gases,  electrochemical  reactions  and  the  transport  of  liquid 
water  through  the  membrane.  Large  scale  and  fine  mesh  sim¬ 
ulations  for  fuel  cells  remain  relatively  scarce  due  to  the  com¬ 
puter  intensive  nature  of  the  numerical  models.  A  key  feature 
of  the  model  implementation  presented  here  is  parallel  com¬ 
puting  which  allows  simulations  with  finer  mesh  resolution 
and/or  a  large  number  of  computational  nodes. 


2.  Model  description  and  field  equations 

The  model  presented  here  is  a  full  three-dimensional 
model  that  resolves  coupled  transport  processes  in  the  mem¬ 
brane,  catalyst  layer,  gas  diffusion  electrodes  and  reactant 
flow  channels  of  a  PEM  fuel  cell.  The  model  was  imple¬ 
mented  via  a  set  of  user  defined  subroutines  in  a  commer¬ 
cial  CFD  code,  Fluent  6.1  [12].  The  implementation  allows 
simulations  using  parallel  processing.  The  assumptions  and 
governing  equations  are  presented  in  this  section. 

2.7.  Assumptions 

Under  constant  load  conditions,  a  fuel  cell  is  assumed  to 
operates  in  steady  state.  Since  the  gas  streams  in  the  flow 
channels  are  humidified,  hydrogen  and  air  at  low  velocities 
(Reynolds  number),  laminar  flow  and  ideal  gas  behaviour  are 
assumed.  The  complex  nature  of  the  transport  processes  in 
PEM  fuel  cells  precludes  systematic  modelling  of  all  pro¬ 
cesses  in  a  three-dimensional  model,  and  phenomena  that 
are  second  order  under  normal  operating  conditions  are  ne¬ 
glected: 

•  All  water  produced  in  the  electrochemical  reactions  is  as¬ 
sumed  to  be  in  the  gas  phase,  and  phase  change  and  two 
phase-transport  are  not  considered. 

•  Dilute  solution  theory  is  used  to  determine  the  species 
diffusion. 

•  The  membrane  is  assumed  to  be  fully  humidified  and  its 
protonic  conductivity  is  taken  to  be  constant. 

•  The  membrane  is  considered  impermeable  to  gases  and 
cross-over  of  reactant  gases  is  neglected  [13]. 


•  Ohmic  heating  in  the  bipolar  plates  and  in  the  gas  diffusion 
electrodes  is  neglected  due  to  high  conductivity. 

•  Ohmic  heating  is  neglected  in  the  membrane.  Heat  trans¬ 
port  in  the  membrane  is  assumed  to  take  place  due  to  con¬ 
duction  only. 

•  Electro-neutrality  prevails  inside  the  membrane.  The  pro¬ 
ton  concentration  in  the  membrane  is  assumed  to  be  con¬ 
stant  and  equal  to  the  concentration  of  fixed  sulfonic  acid 
groups. 

•  The  overpotential  at  the  anode  is  assumed  to  be  constant. 

•  The  gas  diffusion  layer  is  assumed  to  be  homogeneous  and 
isotropic. 


2.2.  Model  equations  in  gas  flow  channels 


CFD  codes  are  structured  around  numerical  algorithms 
that  solve  equations  representing  physical  conservation  laws 
in  order  to  obtain  the  velocity  field  and  associated  mass,  heat 
and  scalar  transport  quantities.  The  three-dimensional  steady 
state  mass  conservation  (continuity)  equation  is  given  by: 


d(pu)  d(pv )  djpw)  _  o 


dx 


dy 


3  z 


where  p  is  the  density  of  the  mixture,  and  u ,  v  and  w  are  the 
velocities  in  the  x,  y  and  z  direction,  respectively.  The  density 
of  the  mixture  is  calculated  using: 


1 

£,■  yt/Pi 


where  is  the  mass  fraction  of  species  i.  The  density  of  each 
species,  pi  is  obtained  from  the  perfect  gas  law  relation: 


PopMi 

RT 


in  which  pop  corresponds  to  the  anode  or  cathode  side  pres¬ 
sure,  Mi  is  the  molecular  weight,  T  is  the  temperature  and  R 
is  the  universal  gas  constant. 

The  flow  field  is  governed  by  the  steady  state  Navier- 
Stokes  equations,  which  express  momentum  conservation  for 
a  Newtonian  fluid.  A  convenient  formulation  for  CFD  is: 


V  •  (pvv)  =  -Vp  +  Vr  +  Sm om  (4) 

where  v  and  p  are  the  velocity  and  pressure  vectors,  r  is  the 
viscous  stress  tensor  matrix  and  Smom  is  a  momentum  source. 


2.2.7.  Mass  transport  equations 

The  steady  state  mass  transport  equation  can  be  written  in 
the  following  general  form: 

V  •  (pvyO  =  —V  •  j;  +  Si  (5) 

where  yi  is  the  mass  fraction  of  species  /,  j,  is  the  diffusive 
mass  flux  vector.  The  source  term  Si  is  set  to  zero  in  the  flow 
channels,  where  no  reaction  takes  place.  Multicomponent  ef¬ 
fects  are  usually  considered  to  be  small  in  dilute  solutions 
and  in  solutions  where  the  species  are  of  similar  size  and 
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nature  [14,15].  In  order  to  account  for  the  tortuosity  in  the 
porous  media  in  the  electrodes,  the  dilute  diffusion  model  had 
to  be  used.  In  Fluent  6.1  there  is  no  straight  forward  way  of 
taking  into  account  the  tortuosity  when  the  Maxwell-Stefan 
equations  are  applied. 

The  diffusive  mass  flux  vector  j*  can  be  written  as  [16]: 


N-l 

it  =  -  PDiivyj 

7=1 


where  Dq  is  the  binary-diffusion  coefficient. 

The  binary  diffusion  coefficients  are  dependent  on  tem¬ 
perature  and  pressure.  They  can  be  calculated  according  to 
the  empirical  relation  [14]: 


Dij  = 


TL15(l/Mi  +  1  /Mf!2 
p  ((E/c  vki)ip  +  (E/t  vkfp) 


x  10 


-3 


where  Dq  is  the  binary  diffusion  coefficient,  T  is  the  temper¬ 
ature  in  Kelvin,  p  is  the  pressure  in  atm,  Mi  is  the  molecular 
weight  of  species  i,  and  Vu  is  the  atomic  diffusion  volume. 
The  values  for  Vki  is  given  by  Cussler  [14]. 


2.2.2.  Energy  transport 

The  energy  equation,  which  expresses  the  first  law  of  ther¬ 
modynamics,  can  be  expressed  as: 


V  •  (y(pE  +  p))  =  V 


fceffVr  -  Ej  h5j  +  (reff  '  v)  )  +v 


where  E  is  the  total  energy,  keff  is  the  effective  conductivity, 
and  j  j  is  the  diffusion  flux  of  species  j  and  h  is  the  enthalpy, 
reff  is  the  effective  stress  tensor  matrix  and  Sh  is  the  source 
term  per  unit  volume  per  unit  time.  However,  the  dissipation 
energy  will  be  very  low  in  the  fuel  cell  due  to  low  velocity 
laminar  flow  and  can  be  omitted  from  the  energy  equation. 
The  first  three  terms  on  the  right-hand  side  of  Eq.  (8)  represent 
the  energy  transfer  due  to  conduction,  species  diffusion,  and 
viscous  dissipation,  respectively. 


2.3.  Model  equations  in  gas  diffusion  layers 


The  gas  diffusion  electrodes  consist  of  carbon  cloth  or 
carbon  fibre  paper  and  can  be  considered  as  porous  media 
through  which  reactant  gases  are  distributed  to  the  catalyst 
layer,  while  the  solid  matrix  of  the  layer  collects  current  and 
connects  the  reaction  zone  to  the  collector  plates.  Fig.  1  shows 
a  schematic  drawing  of  the  cathode  side  of  the  fuel  cell  model, 
illustrating  the  transport  of  reactants,  electrons  and  products. 

The  equations  that  govern  gas  transport  phenomena  in  the 
diffusion  layers  are  similar  to  those  used  in  the  channels  with 
the  addition  of  a  porosity  parameter.  The  mass  conservation 
equation  is  expressed  as: 

(9) 
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Fig.  1.  Schematic  of  the  cathode  side. 

in  which  y  is  the  porosity  and  Sm  is  a  mass  source,  specified 
for  the  anode  and  cathode  according  to  Eqs.  (16)  and  (17) 
presented  in  Section  2.4. 

In  Fluent  6. 1 ,  a  superficial  velocity,  based  on  the  volumet¬ 
ric  flow  rate,  is  used  inside  the  porous  medium.  This  superfi¬ 
cial  velocity  is  also  used  to  ensure  continuity  of  the  velocity 
vectors  across  the  porous  medium  interface.  More  accurate 
simulations  of  porous  media  flows  would  require  solution  of 
the  physical  velocity  throughout  the  flow-field,  rather  than 
the  superficial  velocity.  Fluent  6.1  provides  the  possibility 
of  solving  the  transport  equation  in  the  porous  media  using 
the  physical  velocity.  However,  the  inlet  mass  flow  is  calcu¬ 
lated  from  the  superficial  velocity  and  therefore  the  pressure 
drop  across  the  porous  media  would  be  the  same  whether  the 
physical  or  superficial  velocity  formulation  is  used. 

The  porous  media  model  essentially  consists  of  an  extra 
momentum  sink  term  added  to  the  standard  fluid  flow  Eq. 
(4).  This  source  term  consist  of  two  loss  terms,  viscous  and 
inertial: 

5mom, i  =  —  j  Dijdvj  +  Cij-pVjVj  J  (10) 

V=1  j=l  ) 

where  Smomj  is  the  source  term  for  the  momentum  equation 
in  the  x,  y  or  z  direction,  /x  is  the  dynamic  viscosity,  C  is  the 
inertial  resistance  factor  matrix  and  D  is  a  matrix  containing 
the  inverse  of  the  permeability,  a.  Note  that  the  inertial  loss 
term  is  only  significant  for  high  flow  velocities,  and  its  con¬ 
tribution  is  negligible  in  the  present  model.  The  momentum 
sink  generates  a  pressure  gradient  in  the  porous  region.  The 
steady  state  momentum  equation  in  the  porous  media  using 
the  physical  velocity  formulation  can  then  be  written  as: 

V  •  (ypyy)  =  -yV p  +  Vyr  -  — v  +  Smom  (1 1) 

a 

where  y  and  a  are  porosity  and  permeability  respectively,  and 
the  other  variables  are  defined  as  per  Eq.  (4). 


2.3.1.  Mass  transport  equation  in  porous  media 

The  steady  state  species  transport  equation  in  porous  me¬ 
dia  takes  the  form: 

v  •  ( pyyyi )  =  v  •  {pyDfVyi)  +  St 


V  •  (pyy)  =  Sm 


(12) 
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where  the  source  term  is  given  by  Eqs.  (16)  and  (17)  for  the 
anode  and  cathode  gas  diffusion  layers,  respectively.  Dff  is 
an  effective  diffusion  coefficient  that  takes  into  account  the 
effect  of  additional  drag  by  the  irregular  shape  and  the  actual 
length  of  the  pores  in  comparison  with  a  bundle  of  straight 
parallel  capillaries  with  constant  diameter,  and  is  given  by 
[14,17]: 


(13) 


where  Dt  is  the  diffusion  coefficient,  y  is  the  gas-phase  poros¬ 
ity  and  /ip  is  the  tortuosity  factor.  The  tortuosity  factor  can 
be  split  into  a  factor  that  accounts  for  the  actual  length  of  the 
channel  /il,  and  a  shape  factor  pip ,  (dp  =  PLP2p)  [17].  It 
should  be  noted  that  the  alternative  Bruggemann  correction 
which  is  frequently  used  is  essentially  equivalent  for  low  tor¬ 
tuosities  and  porosities  in  the  range  0.4-0. 5.  In  this  work,  a 
tortuosity  factor  of  3  is  used  [14]. 


2.3.2.  Energy  transport  in  porous  media 

Eq.  (8)  is  used  to  calculate  the  energy  transport  in  the 
porous  media  with  an  effective  thermal  conductivity,  keff> 
calculated  as  the  volume  average  of  the  fluid  conductivity 
and  the  solid  conductivity  (assuming  thermodynamic  equi¬ 
librium),  i.e.: 

&eff  =  ykf  +  ( 1  -  y)ks  (14) 

where  kf  is  the  thermal  conductivity  in  the  fluid  phase  and  ks 
is  the  solid  medium  thermal  conductivity. 

2.3.3.  Potential 

The  potential  distribution  in  the  gas  diffusion  layer  can  be 
calculated  by  applying  the  generic  transport  equation  without 
the  convective  terms. 

-V  •  (orV0)  =  (15) 

where  a  is  the  electronic  conductivity.  The  source  term,  5^, 
is  equal  to  the  local  current  production  which  is  given  in  Eq. 
(19),  and  <fi  is  the  local  potential. 


source  terms  need  to  be  included  in  both  mass  conservation 
and  the  transport  equations. 


2.4.2.  Current  calculations 

The  current  density  at  both  cathode  and  anode  is  calculated 
using  the  Butler- Volmer  equation  [18]: 

•  .-U  .  —  • 

i  =  r  +  i  =  i() 

(19) 


exp 


/3nFr]\ 


RT 


■ 


exp 


(1  —  P)nFrj\ 


RT 


■ 


where  z'o  is  the  exchange  current  density,  n  is  the  number  of 
electrons  per  mole  of  reactant,  rj  is  the  local  overpotential  and 
R  is  the  universal  gas  constant.  /3  is  the  asymmetry  parameter, 
which  is  determined  empirically  to  be  between  0.4  and  0.6 
[18].  The  sensitivity  of  the  model  to  this  parameter  will  be 
discussed  subsequently. 

The  exchange  current  density  is  calculated  from  the  fol¬ 
lowing  relations  [13]: 

Cathode: 


.+  .ref+ 
*0  =  *0 


ro2 


Anode: 


(20) 


(21) 


where  c;  is  the  concentration  of  species  i,  the  superscript  “ref” 
indicates  a  reference  state,  and  y  is  an  empirically  determined 
concentration  parameter.  For  the  cathode  yo2  =  1/2  and 
yH+  =  1/2  [13];  and  for  the  anode  yn2  =  1/4  and  yH+  =  2 
[13]. 

In  this  model  the  concentration  of  protons  at  the  catalyst 
surface  of  the  cathode  and  anode  is  assumed  constant.  Eqs. 
(20)  and  (21)  can  then  be  simplified: 


f  =  kc(co2)YOHcH+)ru+  (22) 

io  =  UcH2)m*  (23) 


2.4.  Catalyst  layer 


2.4.1.  Mass  sources  and  sinks 

The  local  volumetric  source  and  sink  terms  associated  with 
the  electrochemical  reactions  are  proportional  to  the  local 
current  density: 


Sh2  (kgs-1  m-3)  =  - 
So2  (kgs-1  m-3)  =  - 
Sh2o  (kg  s-1  m-3)  = 


M\\2  . 

- -1 

2  E 

(16) 

Mo2  . 

- -i 

4  E 

(17) 

Mu2o  . 

2  E  1 

(18) 

where  Mj  are  the  molar  masses  to  be,  and  F  is  Faraday’s 
constant.  Product  water  is  assumed  in  vapour  form.  These 


where  kc  and  ka  are  experimentally  determined  constants  that 
are  dependent  on  the  geometry  of  the  catalyst  layer. 

No  additional  terms  are  required  for  mass  transport  losses, 
since  these  are  accounted  for  in  the  exchange  current  density 
through  the  concentration  term.  When  the  overpotential  due 
to  mass  transport  losses  is  high,  the  concentration  of  reactants 
in  the  catalyst  layer  is  low,  and  this  results  in  a  low  local 
current  density. 

The  dependence  of  the  exchange  current  density  on  tem¬ 
perature  is  accounted  for  using  the  relation  of  Parthasarathy 
et  al.  (quoted  by  Fuller  and  Newman  [19]): 


io  =  *o(7ref)  exp 


A  E 
~R 


(24) 


where  AE  is  the  activation  energy  (A E  =  73.2 kJ mol  l), 
and  Tref  a  reference  temperature. 
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The  activation  overpotential  stems  from  losses  that  are  as¬ 
sociated  with  the  kinetics  of  the  reactions  taking  place  on  the 
surface  of  the  electrode.  Activation  losses  are  the  most  impor¬ 
tant  irreversibility  in  low  and  medium  temperature  fuel  cells 
[20] .  When  the  reaction  is  not  sufficiently  fast  to  maintain  the 
system  in  equilibrium  at  the  surface,  an  activation  overpoten¬ 
tial  needs  to  be  taken  into  account.  In  many  previous  studies, 
a  constant  surface  overpotential  is  assumed  in  order  to  lin¬ 
earize  the  Butler- Volmer  equation.  In  this  work,  the  local 
activation  overpotential  at  the  cathode  is  obtained  from: 


The  total  cell  potential  for  the  fuel  cell  is  calculated  in  the 
post-processing: 

E  —  ^cathode  ^anode  ^mem  ^contact  (28) 

where  £Cathode  is  set  equal  to  Vinput,cat>  the  input  cathode 
half-cell  potential,  £anode  is  the  anodic  half-cell  potential 
(reversible  cell  potential  minus  overpotentials  at  the  anode), 
77mem  is  the  potential  loss  in  the  membrane,  and  Contact  is  the 
contact  resistance. 


7? act  —  Er  ir  Tinput,cat  (25) 

where  ir  represents  the  ohmic  losses.  Vinput,cat  is  the  half-cell 
potential  at  the  cathode  side,  which  is  an  input  parameter 
used  in  the  model  to  specify  the  operating  point  of  the  cell. 
The  reversible  cell  potential  in  Eq.  (25),  Er,  is  obtained  from 
the  Nernst  equation: 

(26) 

«h*  x / 

where  E°  is  the  open  circuit  voltage  at  standard  pressure  and 
a  is  the  activity.  It  is  assumed  that  the  electrons  are  in  their 
standard  state  and  hence  the  activity  is  taken  as  1 . 

The  ohmic  losses  in  Eq.  (25)  are  obtained  from  solving 
the  potential  equation  in  the  gas  diffusion  layer. 

The  activity  of  the  protons  in  the  membrane  is  difficult  to 
specify.  The  concentration  of  protons  is  assumed  to  be  con¬ 
stant  in  the  membrane.  This  implies  that  the  activity  of  the  two 
half-cell  reactions  should  be  equal,  and  that  the  value  of  the 
proton  activity  will  only  influence  the  local  heat  production 
in  the  half-cell  reactions.  Bernardi  and  Verbrugge  [13]  refer 
to  experimental  measurements  where  the  fixed- charge- site 
concentration  in  the  membrane  is  found  to  be  1 .2  mol  dm-3, 
which  is  almost  equal  to  a  1  molar  solution.  In  Ref.  [21]  the 
activity  of  protons  in  a  1  molar  sulfuric  acid,  H2SO4,  solu¬ 
tion  at  25  °C  is  given  as  0.1316.  Compared  to  other  1  molar 
solutions  this  is  a  low  value  and  it  is  considered  to  be  a  worst 
case  scenario.  In  this  model  the  activity  of  protons  is  assumed 
constant  and  independent  of  temperature. 

The  activity  of  water  vapour  is  given  by: 


2.4.3.  Heat  sources  and  sinks 

The  heat  generated  in  a  fuel  cell  is  due  to  changes  of  en¬ 
thalpy  and  irreversibilities  related  to  charge  transfer.  In  or¬ 
der  to  calculate  the  heat  generation,  both  charged  and  non- 
charged  species  need  to  be  considered.  Values  for  the  entropy 
of  electrons  and  protons  are  required  for  the  computation. 
Lampinen  and  Fomino  [22,23]  presented  a  method  of  calcu¬ 
lating  AG,  AH  and  AS  for  each  half-cell  reaction.  Imple¬ 
menting  this  approach  allows  resolution  of  the  heat  gener¬ 
ation  between  the  cathode  and  anode,  and  a  more  accurate 
prediction  of  the  temperature  profile  in  the  fuel  cell. 

The  energy  balance  for  a  half-cell  can  be  written  as  [22] : 

Q  =  rAH  +  Eel  (29) 

where  Q  is  the  heat  absorbed,  r  is  the  reaction  rate  of  the  half¬ 
cell  reaction,  AH  is  the  half-cell  reaction  enthalpy  and  PQ\  is 
the  electric  power.  The  half-cell  heat  production/absorption 
for  a  real  process  can  be  shown  to  be  [22,24]: 

q=  \i/nF\(AH  +  (—AG))  —  \i\\rj\  =  \i/nF\(TAS)~  \i\fi\ 

(30) 

where  n  is  the  number  of  electrons  per  mole  of  reactants, 
AH,  AG  and  AS  are  the  half-cell  changes  in  enthalpy,  Gibbs 
free  energy  and  entropy,  respectively,  i  is  the  local  current 
density,  r)  is  the  voltage  drop  (overpotential)  due  to  ohmic 
losses  and  reaction  resistance,  and  F  is  Faraday’s  constant. 
The  entropy  change  at  standard  state  with  platinum  catalyst 
is  taken  as  AS  =  0.104  J  mol-1  K-1  for  the  anode  side,  and 
AS  =  —326.36  Jmol-1  K-1  for  the  cathode  side  [22,24]. 


Er  =  E°  - 


RT 

nF 


In 


a 


0.5 

O2 


X 


au2o  = 


P  h2o 

i?w,sat 


(27) 


where  pu2o  is  the  partial  pressure  of  water  vapour  and  jow,sat 
is  the  saturated  water  pressure. 

The  anodic  activation  overpotential  is  much  smaller  than 
the  overpotential  at  the  cathode  [20].  Its  local  variation  is 
therefore  neglected,  and  a  constant  anodic  overpotential  is 
assumed  in  the  present  model.  Due  to  current  conservation, 
the  average  current  density  at  the  anode  and  cathode  have 
to  be  equal.  Hence  the  anodic  overpotential  can  be  found  by 
applying  an  algorithm  which  calculates  the  Butler- Volmer 
equation  using  a  value  for  the  overpotential  that  provides 
an  average  current  density  at  the  anode  equal  to  that  at  the 
cathode. 


2.5.  Membrane 

The  polymer  electrolyte  membrane  acts  primarily  as  a  ma¬ 
terial  barrier  between  the  anode  and  cathode  reactants,  and  as 
a  protonic  conductor.  Both  water  and  heat  are  also  transported 
through  the  membrane. 

2.5.1.  Potential  field  and  heat  transfer 

Electro-neutrality  is  assumed  to  prevail  inside  the  mem¬ 
brane.  Following  Bernardi  and  Verbrugge  [25],  proton  con¬ 
centration  within  the  membrane  is  assumed  to  be  constant 
and  equal  to  the  concentration  of  fixed  sulfonic  acid  groups. 
This  is  a  valid  assumption  for  a  fully  humidified  membrane, 
however,  in  future  work,  approaches  similar  to  those  used  for 
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calculating  charges  in  semiconductors  should  be  considered 
for  computing  the  concentration  of  protons  in  the  membrane. 

Heat  conduction  through  the  membrane  is  represented  by 
Eq.  (8),  without  any  convective  or  source  terms. 

The  potential  field  in  the  membrane  is  computed  using 


-V  •  (/cV0m)  =  0 


(31) 


where  k  is  the  ionic  conductivity,  and  0m  is  the  ionic  potential 
in  the  membrane.  Eq.  (31)  is  solved  subject  to  a  constant 
potential  at  the  anode  and  a  distributed  flux  at  the  cathode 
(see  Section  3.3  for  description  of  boundary  conditions). 

The  ionic  conductivity  is  dependant  on  water  content,  A, 
in  the  membrane.  For  a  fully  humidified  membrane  values 
between  14  and  16.8  are  reported  [26].  In  this  model  a  con¬ 
stant  value  of  15  is  used,  and  the  conductivity  is  obtained 
using  the  empirical  model  of  Springer  et  al.  [26],  yielding  a 
conductivity  of  13.375  S  m-1. 

2.5.2.  Water  transport  in  the  membrane 

The  current  model  considers  only  operation  under  fully 
humidified  conditions.  It  is  common  to  distinguish  between 
three  water  transport  mechanisms  in  polymer  membranes: 
electro-osmosis,  diffusion  and  convection.  The  model  used 
here  is  based  on  the  approach  of  Janssen  [28],  where  rather 
than  prescribing  specific  transport  mechanisms  individually, 
fundamental  thermodynamic  considerations  are  used  to  ac¬ 
count  for  all  the  different  mechanisms  through  a  chemical 
potential  [28]: 


j  —  E+wV/x 


w 


Nw  —  Ew+V  (pm  EWWV  fi 


w 


(32) 

(33) 


where  j  is  the  current  density  of  protons,  0m  is  the  local  po¬ 
tential,  /xw  is  the  local  chemical  potential  of  water  and  the  L 
variables  are  Onsager  coefficients. 

The  physical  interpretation  of  L++  is  the  specific  pro¬ 
ton  conductivity,  k  [28].  According  to  Onsager  relations 
L+w  =  Lw+.  It  can  be  shown  that  Lw+  =  ^  [28],  in  which  § 
is  the  number  of  water  molecules  transported  with  each  pro¬ 
ton.  Further  discussion  and  physical  interpretation  of  Lww  is 
provided  in  Sivertsen  [24]. 

In  Janssen’s  model,  Eqs.  (32)  and  (33)  are  combined  in 
order  to  express  the  water  flux  in  terms  of  the  current  density 
of  protons: 


Nw  —  twj  6ncm  ^  M 


w 


(34) 


in  which 


t  w  — 
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w+ 
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(35) 


(36) 


++ 


The  parameter  tw  is  a  measure  of  the  electro-osmotic  drag  and 
the  term  /memV/xw  is  related  to  the  back  transport  of  water 
in  the  membrane.  Values  for  /mem  for  different  membranes 


are  given  by  Janssen  [28].  In  this  model  /mem  is  assumed 
to  be  constant.  This  is  consistent  with  the  fully  humidified 
membrane  assumption  [28]. 


3.  Boundary  conditions 

Boundary  conditions  are  required  at  all  boundaries  of  the 
computational  domains,  as  well  as  at  internal  interfaces. 

3.1.  Inlet  conditions 

At  the  inlet  of  both  anode  and  cathode  flow  channels  the 
boundary  values  are  prescribed  from  the  stoichiometric  flow 
rate,  temperature  and  mass  fractions.  The  mass  flow  rate 
at  the  inlet  is  prescribed  in  conjunction  with  a  fully  devel¬ 
oped  laminar  flow  profile.  Exact  solutions  for  such  profiles 
are  available  for  a  variety  of  cross-sectional  areas,  but  we 
found  it  computationally  more  effective  to  use  an  approxima¬ 
tion  which  provides  values  within  1%  of  the  exact  solution 
[29]: 
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(38) 


where  wmax  is  the  maximum  velocity,  a  is  the  half-width  of  a 
non-circular  duct,  b  is  the  half-height  of  a  non-circular  duct 
and  um  is  the  average  velocity.  Relations  for  the  values  m  and 
n  are  [29] : 


m  = 


i7+o5© 


2, 


1.4 


n  = 


for  —  <7  J_ 
1U1  b  —  3 


2  +  0.3(f-I),forf  >| 


The  model  was  implemented  to  allow  the  option  of  having 
fully  humidified  inlet  reactant  flows.  The  saturation  pressure 
for  water  vapour  is  given  as  a  quartic  equation  developed  us¬ 
ing  a  regression  of  tabulated  values  from  Moran  and  Shapiro 
[30]. 

pw> sat  (Pa)  =  (1.268366  x  10“8  x  TA  -  1.498267  x  10“5 

xT3  +  0.0067091643  x  T2  -  1.348318703 


xT  +  102.5034101)  x  10' 


(39) 


This  equation  is  valid  for  temperatures  from  323.15  to 
383. 15  K  (50-1 10  °C). 

3.2.  Outlet  conditions 

The  momentum  equation  solver  in  Fluent  uses  a  pressure 
correction  method  which  does  not  allow  specification  of  two 
different  reference  pressures  in  the  computational  domain. 
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Since  the  reactant  gas  flow  channels  are  separate  and  gener¬ 
ally  at  different  pressures,  pressure  boundary  conditions  are 
used  at  the  outlets. 

3.3.  Boundary  conditions  at  internal  interfaces 

3.3.1.  Interface  between  the  electrode  and  the  flow 
channel 

At  the  interface  between  the  electrode  and  the  flow  chan¬ 
nel,  a  user  defined  function  (UDF)  was  implemented  to  en¬ 
force  a  zero  flux  of  electrons.  This  is  done  by  overwriting  the 
value  of  the  potential  in  the  first  cell  in  the  flow-channel  side 
of  the  interface  with  the  same  value  as  the  adjacent  cell  on 
the  electrode  side,  thereby  enforcing  a  zero  potential  gradient 
normal  to  the  interface. 

3.3.2.  Interface  between  the  electrode  and  the  membrane 

Due  to  the  structure  of  the  Fluent  CFD  code,  the  inter¬ 
face  between  the  membrane  and  the  electrodes  is  defined 
as  a  wall  (impermeable  boundary).  This  is  done  mainly  to 
prevent  any  crossover  of  species  and  electrons  through  the 
membrane,  but  also  to  prevent  pressure  related  problems. 
The  wall  has  a  fluid  region  on  each  side,  and  each  side  is 
therefore  treated  as  a  distinct  wall.  This  is  implemented  by 
creating  a  “shadow”  of  the  wall  cell  layer.  The  “walls”  are 
thermally  coupled,  such  that  heat  transfer  can  be  directly 
computed  from  the  solution  in  the  adjacent  cells,  and  hence 
no  additional  boundary  conditions  is  required  for  the  energy 
equation.  The  use  of  an  internal  wall  is  an  expedient  way 
of  implementing  our  model  in  Fluent  and  it  is  physically 
consistent,  but  more  straightforward  alternatives  should  be 
explored. 

The  potential  profile  in  the  membrane  is  solved  by  setting 
the  ionic  potential  to  zero  at  the  anode  side  interface  between 
the  catalyst  and  the  membrane.  This  is  based  on  the  assump¬ 
tion  that  the  potential  at  the  anode  side  is  uniform.  At  the 
membrane-catalyst  interface  on  the  cathode  side,  a  negative 
flux  equal  to  the  consumption  of  protons  in  the  adjacent  cell 
in  the  catalyst  layer  is  specified. 

As  described  earlier,  the  gradient  of  the  chemical  potential 
is  used  to  calculate  the  water  flux  through  the  membrane.  The 
profile  of  the  chemical  potential  in  the  membrane  is  computed 
using  the  water  transport  equation  subject  to  specified  values 
at  the  membrane/electrode  interfaces.  Local  equilibrium  is 
assumed  at  the  interfaces,  i.e.  the  water  in  the  membrane 
and  the  water  vapour  in  the  electrode  are  in  equilibrium.  The 
boundary  conditions  at  the  anode  and  cathode  interfaces  are 
[28]: 

/C”  =  nl  +  RT  In  ^  (40) 

pKJ 

where  /z™em  is  the  chemical  potential  of  water  vapour  in  the 
membrane  at  the  interface,  /z^  is  the  chemical  potential  at 
standard  conditions,  pw  is  the  water  vapour  pressure  in  the 
electrode  at  the  interface,  and  p°  is  the  standard  pressure. 


3.4.  Land  areas 

In  the  areas  where  the  gas  diffusion  electrodes  are  in  con¬ 
tact  with  the  bipolar  plates  (the  land  area)  a  constant  reference 
voltage  equal  to  zero  is  applied  as  a  boundary  condition  both 
at  the  anode  and  at  the  cathode.  At  all  other  walls,  the  electron 
flux  is  set  to  zero. 

3.5.  Walls 

On  all  walls  the  no-slip  boundary  condition  is  applied  for 
the  momentum  equations. 

Different  boundary  conditions  for  the  energy  equation  can 
be  employed:  specified  heat  flux,  temperature  or  convective 
heat  transfer.  The  specific  boundary  condition  will  be  pre¬ 
sented  with  each  of  the  cases  investigated  in  this  paper. 

4.  Computational  procedure 

The  governing  transport  equations  are  solved  subject  to 
the  various  boundary  conditions  presented  in  the  preceding 
section  using  the  Fluent  6.1  [12]  CFD  code.  Fluent  6.1  is 
a  parallel  code  using  a  finite  volume  method  and  an  itera¬ 
tive  segregated  implicit  solver.  Second  order  discretization 
schemes  are  utilized  for  all  transport  equations  (i.e.  central 
differencing  for  diffusion  terms,  and  second  order  upwind  for 
convection  terms). 

The  model  equations  specific  to  physico-chemical  trans¬ 
port  in  PEM  fuel  cells  presented  earlier  were  implemented 
using  a  set  of  UDFs  written  in  C  and  dynamically  linked  with 
the  Fluent  source  code.  These  UDFs  were  developed  such  as 
to  maintain  the  parallel  computing  capabilities  which  sig¬ 
nificantly  enhances  the  useability  of  the  code  for  problems 
necessitating  a  large  number  of  computational  nodes.  An¬ 
other  important  feature  of  the  model  presented  here  is  an 
iterative  voltage-current  algorithm  also  implemented  using 
UDFs.  This  algorithm  allows  prediction  of  the  cell  current 
density  based  on  a  target  cell  operating  voltage,  and  does  not 
require  the  assumption  of  a  specific  constant  overpotential 
over  the  entire  cell,  but  rather  allows  prediction  of  the  spatial 
variation  of  the  local  surface  overpotential. 

5.  Test  case  geometry  and  parameters 

The  geometry,  parameters  and  operating  conditions  used 
in  the  simulations  correspond  to  the  experimental  test  case 
of  Wang  et  al.  [31].  It  is  common  to  compare  results  from 
model  predictions  with  polarization  curves  obtained  from 
experiments  or  other  models.  This  is  not  a  very  satisfactory 
way  of  verifying  model  performance.  Until  progress  is  made 
in  the  difficult  problem  of  obtaining  detailed  and  reliable  in 
situ  experimental  data,  we  will  continue  to  use  polarization 
curves  and  complement  the  assessment  with  an  analysis  of 
the  physical  results  obtained  with  the  model.  Given  the  state 
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of  development  of  PEM  fuel  cell  modelling,  the  complexity 
of  the  transport  phenomena,  and  the  large  range  of  physical 
scales  involved,  and  the  uncertainty  in  determining  some  of 
the  physico-chemical  parameters,  it  is  any  case  unrealistic  to 
expect  accurate  quantitative  agreement  between  experiments 
and  computational  model  simulations.  Simulation  based  on 
physically  representative  models  should  however  yield  cor¬ 
rect  relative  trends  and  provide  valuable  insight  and  guidance 
for  design  and  optimization. 

The  geometry  of  the  fuel  cell  simulated  in  this  work  is  a 
straight  section  consisting  of  bipolar  plates  with  flow  chan¬ 
nels  separated  by  a  membrane-electrode  assembly.  The  di¬ 
mensions  of  the  fuel  cell  section  are  given  in  Table  1  and 
correspond  to  the  experimental  test  case  of  Ref.  [31].  A  non- 
uniform  grid  was  used  to  minimize  the  computational  re¬ 
quirements  while  allowing  proper  resolution  in  high  gradient 
regions  (near  wall  regions;  porous  electrodes).  It  should  be 
noted  that  though  the  implementation  allows  for  the  resolu¬ 
tion  of  the  catalyst  layer,  the  simulations  presented  here  use 
a  single  computational  layer  for  the  catalyst  layer. 

A  grid  sensitivity  study  was  performed  to  establish  the 
grid  resolution  requirements.  Three  meshes  with  546,000 
(grid  A),  155,925  (grid  B)  and  66,240  (grid  C)  computational 
cells  each  were  generated.  At  a  cell  potential  of  0.798  V,  the 
coarsest  grid  C  and  the  medium  grid  B  yield  deviations  of 
about  1.2%  and  0.9%  for  the  average  current  density  com¬ 
pared  to  the  fine  grid  A.  All  simulations  presented  here  were 
performed  using  the  finest  grid  in  which  the  anode  and  the 
cathode  flow  channel  are  each  divided  into  20  x  18  x  150 
grid  cells;  the  membrane  into  40  x  7  x  150  computational 
cells;  and  the  gas  diffusion  layers  into  40  x  12  x  150  cells. 
The  catalyst  layers  thickness  consists  of  one  computational 
volume,  i.e.  a  total  of  40  x  1  x  150  cell  each.  It  should  be 
noted  that  the  parallel  solver  allows  us  to  use  a  large  compu¬ 
tational  grid  (total  of  546,000  cells)  that  allows  much  finer 
resolution  than  those  employed  in  previous  work  (e.g.  [5,6]). 

Table  2  provides  the  basic  experimental  operating  condi¬ 
tions  [31]  used  for  the  simulations.  Since  the  computational 


Table  1 


Physical  dimensions  for  the  straight  channel  fuel  cell  section 

Parameter 

Value 

Unit 

Cell  width 

2.0  x  10“3 

m 

Channel  length 

0.07 

m 

Channel  height 

1.0  X  10“3 

m 

Channel  width 

1.0  x  10“3 

m 

Land  area  width 

1.0  x  10“3 

m 

Electrode  thickness 

0.3  x  10“3 

m 

Catalyst  layer  thickness 

1.29  x  10“5 

m 

Membrane  thickness 

0.108  x  10“3 

m 

Table  2 

Operation  parameters  for  the  straight  channel  fuel  cell  section 

Parameter 

Value 

Unit 

Inlet  temperature,  anode  and  cathode 

80 

°C 

Anode  side  pressure 

3 

atm 

Cathode  side  pressure 

3 

atm 

Anode  stoichiometric  flow  rate 

3 

— 

Cathode  stoichiometric  flow  rate 

3 

— 

Relative  humidity  of  inlet  gases 

100 

% 

Oxygen/nitrogen  ratio 

0.79/0.21 

— 

Table  3 

Inlet  mass  fraction  (fully  humidified  flow) 

Species 

Mass  fraction  (%) 

H2O,  anode 

62.253 

H2,  anode 

37.747 

H2O,  cathode 

10.344 

O2,  cathode 

20.885 

N2,  cathode 

68.771 

domain  does  not  include  a  serpentine  flow  channel,  the  mass 
flow  rate  is  adjusted  to  maintain  a  stoichiometry  of  3.  The  re¬ 
actants  are  assumed  to  be  fully  humidified  at  the  inlets  (Table 

3). 

The  gas  diffusion  layer  and  membrane  properties  used  in 
the  base  case  simulations  are  given  in  (Tables  4  and  5).  In 
the  literature,  the  asymmetry  parameter  /3  (also  referred  to 
as  the  transfer  coefficient)  is  ascribed  values  ranging  from 


Table  4 


Electrode  properties 


Parameter 

Symbols 

Value 

Unit 

Reference 

Electrode  porosity 

y 

0.4 

— 

[31] 

Permeability 

a 

1.76  x  10"11 

2 

in 

[31] 

Tortuosity 

P  p 

3.0 

— 

[14] 

Thermal  conductivity 

Solid  region 

k 

150.6 

W  m_l  K-1 

[6] 

Electronic  conductivity 

0 

100 

Sm'1 

[36] 

Asymmetry  parameter 

Anode 

0.5 

— 

[18] 

Cathode 

1.0/0. 5 

— 

[13,18,32] 

Concentration  parameter 

Hydrogen 

m2 

0.25 

— 

[13] 

Oxygen 

Yo2 

0.5 

— 

[13] 

Constant,  anode 

ka 

17.00  x  107 

— 

— 

Constant,  cathode  f  =  1.0 

kc 

40.75  x  10“n 

— 

— 

Constant,  cathode  f>  =  0.5 

kc 

0.8093958 

— 

— 
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Table  5 


Membrane  properties 


Parameter 

Symbols 

Value 

Unit 

References 

Thermal  conductivity 

k 

0.67 

Win-1  K-1 

[6] 

Proton  conductivity 

K 

13.375 

S  m“ 1 

[26] 

Permeability  in  the  membrane 

^mem 

2.35  x  lO"7 

mol2  sm~3  kg-1 

[28] 

Fixed-charge  concentration 

Number  of  water  molecules 

cH+ 

1200 

mol  m_  1 

[13] 

Transported  per  proton 

$ 

3 

— 

[28,27] 

1  to  0  [32].  Because  of  this  variability,  a  parametric  study 
was  undertaken  to  assess  the  effect  of  the  asymmetry  param¬ 
eter  at  the  cathode,  where  this  parameter  was  varied  for  0.5 
to  1.0. 

It  should  be  noted  that  the  value  given  by  Bernardi  and  Ver- 
brugge  [13]  for  the  transfer  coefficients  at  the  anode,  is  not 
strictly  consistent  with  the  Butler- Volmer  equation  [18,32] 
which  requires  that  the  coefficients  for  the  forward  and  back¬ 
ward  components  of  should  add  to  unity.  The  values  for  trans¬ 
fer  coefficients  at  the  anode  given  by  Bernardi  and  Verbrugge 
[13],  which  also  was  used  by  Wang  et  al.  [31],  is  therefore 
not  used  in  this  model.  However,  since  the  activation  overpo¬ 
tential  at  the  anode  is  much  smaller  than  that  at  the  cathode, 
the  value  of  the  charge  transfer  coefficients  at  the  anode  is  of 
less  importance. 

A  constant  temperature  of  80  °C  is  applied  as  a  boundary 
condition  to  the  outer  walls  of  the  fuel  cell.  The  thermal  con¬ 
ductivity  of  the  bipolar  plates  is  set  to  35  W  m-1  K-1  and  the 
electrical  resistance  is  set  to  50  x  10-6  £2  m,  which  are  a  mid¬ 
range  values  for  the  data  presented  by  Middelman  et  al.  [33]. 
In  the  experiments  of  Wang  et  al.  [31]  gold-plated  copper 
plate  were  used  to  connect  the  bipolar-plates  with  the  exter¬ 
nal  circuit;  the  electrical  contact  resistance  is  consequently 
set  to  zero  in  the  computational  model. 


All  other  properties  used  for  the  fluids  and  solids  are  taken 
from  Mills  [34] . 

6.  Results  and  discussion 

The  results  obtained  for  the  straight  channel  case  is  pre¬ 
sented  and  discussed  in  this  section.  All  computations  were 
preformed  on  a  Linux  computer  cluster  with  dual  2000+ 
AMD  Athlon  processors  on  each  node.  The  maximum  num¬ 
ber  of  processors  used  were  eight.  The  number  of  iterations 
required  and  the  computation  time  was  dependent  on  the  ini¬ 
tial  conditions  specified  and  the  input  of  the  cathode  half-cell 
voltage. 

6.1.  Polarisation  curve  and  effect  of  asymmetry 
parameter 

Good  agreement  between  measured  and  computed  polari¬ 
sation  curves  is  not  sufficient  to  assess  the  predictive  capabil¬ 
ities  of  a  model,  but  is  a  prerequisite.  In  order  to  investigate 
the  effect  of  the  asymmetry  parameter  on  the  polarisation 
curve,  simulations  were  performed  for  (3  =  1.0  and  0.5.  Fig. 
2  compares  the  computed  and  measured  polarisation  curves. 


Fig.  2.  Polarisation  curve:  comparison  of  simulations  and  experiments  [31]. 
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In  the  activation  and  the  ohmic  region  (low  and  mid-range 
current  densities)  of  the  polarisation  curve,  the  results  from 
the  model  are  in  good  agreement  with  experiments,  with  de¬ 
viations  of  less  than  5.0%  for  p  =  1.0  and  less  than  1.1% 
with  p  =  0.5.  At  very  low  current  densities  the  simulations 
with  p  =  0.5  appear  to  follow  the  experimental  results  better 
than  the  simulations  with  p  =  1.0.  This  suggests  the  actual 
value  of  the  asymmetry  parameter  should  be  closer  to  0.5  than 
1 .0,  which  is  consistent  with  the  range  of  values  reported  by 
Hamann  et  al.  [18]. 

As  expected,  the  model  is  unable  to  reproduce  the  experi¬ 
mental  data  at  high  current  densities.  The  discrepancy  in  this 
region  is  attributed  in  part  to  the  assumption  of  single  phase 
transport.  In  practice,  however,  the  formation  of  liquid  is  ex¬ 
pected  to  limit  mass  transport  at  higher  current  densities  [8]. 
Another  effect  that  may  lead  to  a  drop  in  the  cell  voltage  at 
higher  current  densities,  is  the  shift  of  the  reaction  zone  in 
the  catalyst  layer  away  from  the  membrane  interface  at  higher 
currents  [9] .  The  protons  need  to  be  transported  further  out  in 
the  catalyst  layer  due  to  depletion  of  oxygen  in  the  catalyst, 
and  this  lead  to  increased  ohmic  losses.  This  effect  is  not  re¬ 
solved  here  since  the  catalyst  is  modelled  as  a  single  layer  of 
computational  cells. 

Fig.  3  shows  the  local  current  distributions  obtained  in  the 
cathode  catalyst  layer  for  two  values  of  asymmetry  parameter, 
p  =  1.0  and  0.5.  The  p  =  1.0  case  exhibits  markedly  higher 
peaks  in  the  shoulders  near  the  inlet.  The  relative  difference 
between  the  two  case  is  shown  in  Fig.  4.  The  difference  in  the 
current  density  profiles  can  be  explained  by  the  fact  that  for 
any  given  activation  overpotential,  the  Butler- Volmer  equa¬ 
tion  yields  a  higher  current  density  for  a  higher  value  of  the 
asymmetry  parameter,  and  hence  the  current  density  profile 
will  exhibit  larger  maxima.  Since  the  average  current  den¬ 
sity  is  the  same  in  both  cases,  the  total  oxygen  consumption 
is  also  the  same.  Whereas  the  difference  in  the  local  current 
density  is  greater  than  18%  at  some  locations,  the  cathode 
half-cell  potential  and  the  cell  potential  are  virtually  identi¬ 
cal  (maximum  deviation  less  than  0.2%).  This  result  shows 
that  even  though  the  results  obtained  for  the  two  asymmetry 
parameters  yield  identical  data  on  the  polarization  curve,  the 
current  distributions  are  significantly  different.  The  predic¬ 
tion  of  physically  representative  distributions  is  important  in 
practice  as  non-uniformity  can  have  a  significant  impact  on 
the  design  and  longevity  of  a  fuel  cell. 

6.2.  Relative  influence  of  activation  and  ohmic  losses 

Figs.  5  and  6  present  the  current  density,  activation  overpo¬ 
tential  and  oxygen  mass  fraction  distributions  for  two  current 
densities,  /avg  =  0.269  and  0.989  A  cm-2.  For  the  lower  load 
case,  cathode  activation  overpotential  and  the  local  current 
density  are  fairly  uniform,  but  both  vary  much  more  in  the 
higher  load  case.  Both  figures  show  that  the  current  density 
profiles  correlate  with  the  activation  overpotential,  but  not 
with  the  oxygen  concentration.  These  results  are  radically 
different  from  the  intuitively  expected  profiles  and  from  pre- 
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Fig.  3.  Effect  of  asymmetry  parameter  on  cathodic  current  density  distribu¬ 
tion  at  iavg  0.269  A  cm-2.  Top:  ft  =  1.0;  bottom:  f  =  0.5. 


dictions  obtained  with  models  assuming  constant  overpoten¬ 
tial  (e.g.  [6]),  where  the  current  densities  are  highest  in  the 
centre  of  the  channel  and  coincide  with  the  highest  reactant 
concentrations.  The  present  results  are  also  consistent  with 
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Fig.  5.  Isocontours  for  current  density,  activation  overpotential  and  oxy¬ 
gen  mass  fraction  at  the  cathode  catalyst  layer  at  V  =  0.798  V,  iawg  = 
0.269  A  cm-2, 13  =  0.5. 


recent  studies  where  distributed  activation  potential  and  com¬ 
plete  charge  transport  are  accounted  for  [9,35]. 

The  current  density  maxima  under  the  land  areas  is  as¬ 
sociated  with  the  fact  that  ohmic  losses  in  the  gas  diffusion 
layer  influence  the  activity  at  the  catalyst  more  than  the  con¬ 
centration  losses.  The  electric  current  path  from  the  area  of 
the  catalyst  layer  under  the  flow  channel  is  longer  than  the 
path  from  the  area  of  the  catalyst  layer  under  the  land  areas 
(see  Fig.  1).  The  potential  field  in  the  cathodic  and  the  anodic 
gas  diffusion  electrodes  are  shown  in  Fig.  7  at  the  middle 
plane  of  the  cell.  The  isopotential  lines  are  normal  to  the  flow 
channel  and  the  side  walls,  while  there  is  a  gradient  into  the 
land  areas  where  electrons  flow  into  the  bipolar  plate.  The 
distributions  exhibit  gradients  in  both  v  and  z  direction  due 
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gen  mass  fraction  at  the  cathode  catalyst  layer  at  Vref  =  0.652  V,  zavg  = 
0.989  A  cm-2,  yd  =  0.5. 


to  the  non-uniform  local  current  production  and  show  that 
ohmic  losses  are  larger  in  the  area  of  the  catalyst  layer  under 
the  flow  channels.  Note  that  the  potential  along  the  anode  side 
interface  (bottom)  is  assumed  uniform  and  is  the  reference 
potential  (VmQm  =  0  V). 

Since  the  conductivity  of  the  gas  diffusion  layer  is  ex¬ 
pected  to  have  a  strong  influence  on  the  activation  overpo¬ 
tential,  and  hence  on  the  current  distribution,  parametric  sim¬ 
ulations  for  several  conductivities  were  performed  to  assess 
this  effect.  Fig.  8a-c  shows  the  changes  in  the  profiles  when 
the  resistance  in  the  gas  diffusion  layer  is  reduced.  Since  the 
fuel  cell  will  operate  at  different  loads  (fixed  cathode  half-cell 
potential)  when  the  conductivity  in  the  gas  diffusion  elec¬ 
trode  is  changed,  relative  current  density  profiles  are  exam¬ 
ined.  Results  are  presented  for  three  conductivity  values  while 
maintaining  all  other  parameter  constant:  100  Sm-1  [36], 
825  Sm-1  [9]  and  1500Sm_1.  With  increasing  electronic 
conductivity  (or  reduced  electrode  thickness),  ohmic  losses 
through  the  electrode  decrease,  the  concentration  losses  be¬ 
come  larger  relative  to  the  ohmic  losses,  and  the  local  current 
density  maxima  shift  toward  the  centre  of  the  channel  where 
the  losses  are  smallest. 

6.3.  Further  discussion 

The  potential  distribution  in  the  fuel  cell  has  been  shown  to 
govern  the  current  density  distribution.  In  this  study,  effects 
associated  with  two-phase  flow  have  not  been  considered. 
These  might  have  a  large  impact  on  the  actual  current  profile 
at  higher  current  densities,  and  liquid  water  pockets  are  more 
likely  to  form  under  the  land  area  of  the  gas  diffusion  layer 
[37,8].  This  would  inhibit  the  diffusion  of  the  reactant  to 
the  catalyst  layer  and  thereby  reduce  the  current  density  in 
this  area,  thus  counterbalancing  some  of  the  trends  reported 
in  this  study.  For  simulations  corresponding  to  intermediate 
current  densities,  where  condensation  is  not  expected  to  play 
a  significant  role,  the  result  reported  here  are  expected  to  be 
physically  representative. 

In  this  work  the  gas  diffusion  layer  was  assumed  to  be  ho¬ 
mogeneous  and  isotropic,  whereas  actual  carbon  paper  ma¬ 
terials  commonly  used  to  fabricate  gas  diffusion  electrodes 
are  anisotropic,  and  the  permeability  is  likely  to  be  lower 
in  the  plane  of  the  carbon  paper  layer  than  across  it.  The 
species  transport  in  the  in  plane  directions  would  thus  be 
expected  to  be  lower  than  in  the  direction  perpendicular  to 
the  gas  diffusion  layer.  A  potentially  more  important  effect 
of  the  non-isotropic  gas  diffusion  layer  is  that  the  tortuos¬ 
ity  factor  could  be  different  in  the  two  directions,  leading 
to  a  lower  diffusion  into  the  areas  under  the  land  areas.  Ac¬ 
counting  for  the  anisotropic  properties  of  the  gas  diffusion 
layer  would  likely  yield  lower  concentration  of  reactant  gas 
under  the  land  area,  and  hence,  lower  local  current  densi¬ 
ties  in  this  area.  There  are  no  fundamental  difficulties  in  ac¬ 
counting  for  the  anisotropic  properties  in  the  CFD  model,  the 
problem  is  rather  in  obtaining  reliable  data  in  prescribing  the 
properties. 
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Fig.  7.  Cross-section  contours  of  the  potentials  in  the  MEA  at  a  midway 
location  in  the  channel  ( V  =  0.652  V,  iavg  =  0.989  A  cm-2,  ft  =  0.5).  (a) 
Cathode,  (b)  membrane,  and  (c)  anode. 
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Fig.  8.  Plots  of  the  relative  current  densities  (7  local/  zavg)  for  different  conduc¬ 
tivities  Oavg  =  0.269  A  cm-2;  V  =  0.798  f  =  0.5).  Electric  conductivity: 
(a)  100  Sm-1,^)  825  8  m-1,  and  (c)  1500  8  m-1. 


6.4.  Conclusions 

A  three-dimensional  model  of  a  proton  exchange  mem¬ 
brane  fuel  cell  has  been  developed  and  implemented  in  the 
framework  of  a  CFD  code.  The  model  implementation  takes 
advantage  of  the  parallel  processing  architecture  of  the  Fluent 
CFD  code,  thus  allowing  fast  simulations  even  with  fine  grid 
resolution  and/or  large  computational  domains.  The  grid  res¬ 
olution  used  in  this  study  is  significantly  finer  than  those  re¬ 
ported  in  previous  computational  work.  The  model  accounts 
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for  the  fluid  transport  inside  the  channels  and  the  gas  diffusion 
electrodes  as  well  as  heat  transfer.  The  cathodic  overpoten¬ 
tial  distribution  is  resolved,  rather  than  assumed  uniform  and 
constant,  allowing  more  accurate  predictions  of  local  current 
densities. 

Global  comparisons  show  good  agreement  between  the 
model  and  experimental  results.  The  computational  anal¬ 
ysis  shows  that  substantially  different  spatial  distributions 
can  be  obtained  by  varying  the  asymmetry  parameter  with 
no  noticeable  change  in  the  polarization  curves.  This  fur¬ 
ther  highlights  that  global  comparison  between  experimen¬ 
tal  and  predicted  results  is,  as  pointed  out  in  previous 
studies,  insufficient  to  validate  computational  model.  De¬ 
tailed  experimental  data  providing  spatially  resolved  dis¬ 
tributions  of  key  quantities  is  essential  for  proper  valida¬ 
tion. 

The  predicted  distribution  of  current  densities  show  pro¬ 
files  which  are  fundamentally  different  from  the  distribution 
obtained  in  multi-dimensional  models  that  do  not  account 
for  distributed  overpotentials.  The  maximum  current  density 
occurs  under  the  land  areas  as  a  result  of  the  dominant  in¬ 
fluence  of  ohmic  losses  over  concentration  losses  on  the  ac¬ 
tivity  at  the  catalyst  layer.  The  study  shows  that  changing 
the  conductivity  radically  alters  the  current  distribution  by 
changing  the  relative  influence  of  ohmic  to  activation  over¬ 
potentials. 
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